* To save data in MATLAB: save filename.txt -ascii -tabs

*******************************************************************************************************
* Collect the dataset:
*******************************************************************************************************

* Identified shocks from the VAR:

clear
insheet using strtrans_rev.csv, comma
ren v1 shock_SBT
ren v2 shock_n
ren v3 shock_rest
cap drop v*
drop if shock_SBT==.

/*
use varshocks201111
*/

gen nr = _n
gsort -nr
drop nr
local nrobs = _N+9 /* FD + 8 lags */
set obs `nrobs'
gen nr = _n
gsort -nr
drop nr

/*
foreach VAR of varlist shock_GaliOLS* {  /*OLS only 4 lags */
	ren `VAR' temp
	gen `VAR' = temp[_n+4]
	drop temp
}
*/

gen t = _n
drop if t>110 & shock_SBT==.
compress
sort t
save shocks, replace

* Variables used in the VAR (have been seasonally adjusted):

clear
insheet using data_spec26.csv, comma
ren quarter qrt

gen nr = _n
gsort -nr
drop nr
local nrobs = _N+1 /* FD */
set obs `nrobs'
gen nr = _n
gsort -nr
drop nr
replace year = 1979 if year==.
replace qrt = 1 if qrt==.

gen time = year + (qrt-1)/4
sort time
gen t = _n
sort t
merge 1:1 t using shocks
drop _merge
tsset t

* Calculate series in levels
foreach VAR in prem rele prod hours wl {
	gen `VAR' = 0
	replace `VAR' = L.`VAR' + d_`VAR' if L.`VAR'~=.
}

compress
tsset t
save datashocks, replace


*******************************************************************************************************
* Check data:
*******************************************************************************************************

clear
use datashocks

ren prem sp
*rele
*prod
ren hours hrs
ren wl wL

su sp prod hrs wL rele
su D.sp D.prod D.hrs D.wL D.rele

gen wH = wL + sp
gen hL = hrs - 100*ln( 1 + exp(rele/100) )
gen hH = hL + rele

/* OLD DETRENDING METHOD:
gen byte post97 = (year>=1997)
pause off
foreach VAR of varlist sp rele prod hrs wL wH hL hH {
	reg D.`VAR'
	predict d`VAR'_dt_naive, resid
	reg D.`VAR' post97
	predict d`VAR'_dt, resid
	gen `VAR'_dt_naive = 0 in 1
	gen `VAR'_dt = 0 in 1
	replace `VAR'_dt_naive = `VAR'_dt_naive[_n-1]+d`VAR'_dt_naive if d`VAR'_dt_naive~=.
	replace `VAR'_dt = `VAR'_dt[_n-1]+d`VAR'_dt if d`VAR'_dt~=.
	graph7 `VAR'_dt `VAR'_dt_naive time, xlabel ylabel sy(do) c(ll) xline(1997)
	pause
}
drop *_naive
su dsp* drele* D.sp* D.rele*
*/


*******************************************************************************************************
* Identification graphs:
*******************************************************************************************************

* SA, normalize and smooth series and shock for nicer graphs:
*tab qrt, g(q_)
*drop q_1
foreach VAR of varlist sp prod hrs rele {
	cap gen d`VAR' = D.`VAR'
	*reg d`VAR' q_*
	*predict d`VAR'_sa, resid
	*drop d`VAR'
	*ren d`VAR'_sa d`VAR'
	qui su d`VAR'
	replace d`VAR' = (d`VAR' - r(mean))/r(sd)
	gen d`VAR'_sm = (1/4)*L4.d`VAR' + (1/2)*L3.d`VAR' + (3/4)*L2.d`VAR' + L.d`VAR' + d`VAR' + F.d`VAR' + (3/4)*F2.d`VAR' + (1/2)*F3.d`VAR' + (1/4)*F4.d`VAR'
}
foreach VAR of varlist shock* {
	gen `VAR'_sm  = (1/4)*L4.`VAR'  + (1/2)*L3.`VAR'  + (3/4)*L2.`VAR'  + L.`VAR'  + `VAR'  + F.`VAR'  + (3/4)*F2.`VAR'  + (1/2)*F3.`VAR'  + (1/4)*F4.`VAR'
}

/*
* GRAPHS ABOUT GALI IDENTIFICATION:

*graph7 shock_Gali_n shock_n
graph7 dprod_sm dhrs_sm shock_Gali_n_sm time, xlabel ylabel sy(..d) c(ll[.]l)

* Is the extremely close correspondence between dprod and shock_n due to Minnesota prior?
graph7 dprod_sm dhrs_sm shock_GaliOLS_n_sm time, xlabel ylabel sy(..d) c(ll[.]l)
*/

* GRAPHS ABOUT SBT IDENTIFICATION:
graph7 dsp_sm dprod_sm dhrs_sm shock_SBT_sm shock_n_sm time, xlabel ylabel sy(...do) c(ll[-]l[.]ll[-])
graph7 dsp_sm shock_SBT_sm dhrs_sm time, xlabel ylabel sy(...) c(ll[-]l[.])
graph7 dprod_sm shock_n_sm dhrs_sm time, xlabel ylabel sy(...) c(ll[-]l[.])

* Nice graph for letter to the referees:
twoway (line dsp_sm time, lcolor(blue) lwidth(medthin) lpattern(solid)) (line shock_SBT_sm time, lcolor(blue) lwidth(medthin) lpattern(dash)) (line dprod_sm time, lcolor(red) lwidth(medthick) lpattern(solid)) (line shock_n_sm time, lcolor(red) lwidth(medthick) lpattern(dash)) (line dhrs_sm time, lcolor(green) lwidth(medthick) lpattern(dot))
graph save graph_identification, asis replace
graph export graph_identification.eps, as(eps) preview(on) replace
pause

keep t time sp prod hrs rele wL wH hL hH shock*


*******************************************************************************************************
* Check whether Solow residuals similar to identified SBT shocks:
*******************************************************************************************************

local sigma = 1.5
gen dsolow = D.sp + (1/`sigma')*D.rele
*gen dsolow_dt = D.sp_dt + (1/`sigma')*D.rele_dt

qui su dsolow if shock_SBT~=.
replace dsolow = ( dsolow - r(mean) )/r(sd)
*qui su dsolow_dt if shock_SBT~=.
*replace dsolow_dt = ( dsolow_dt - r(mean) )/r(sd)

su dsolow* shock_SBT
su dsolow* shock_SBT if shock_SBT~=.
corr shock_SBT dsolow

forvalues sig = 10 20 to 100 {
	local sigma = `sig'/10
	gen dsolow_`sig' = D.sp + (1/`sigma')*D.rele
	qui su dsolow_`sig' if shock_SBT~=.
	replace dsolow_`sig' = ( dsolow_`sig' - r(mean) )/r(sd)
}
corr shock_SBT dsolow_*

pause on
*graph7 dsolow shock_SBT time, xlabel(1980 1982 to 2006) ylabel sy(do) c(l[.]l)
*pause
*graph7 dsolow shock_SBT, xlabel ylabel 

save datashocks_more, replace
outsheet dsolow shock_SBT using solowres_and_sbtshock.txt, nonames nolabel replace


*******************************************************************************************************
* Construct IRFs as in Basu-Fernald-Kimball:
*******************************************************************************************************

*forvalues maxlag = 2(2)6 {
local maxlag = 6
if 1==1 {

clear
set more off

* Generate dataset for IRFs
set obs 140
gen t = _n-40
gen shock = 0
replace shock=1 if t==0
sort t
tsset t
save irfs, replace
clear

foreach VAR in "sp" "prod" "hrs" "rele" "wL" "wH" "hL" "hH" {
foreach SHOCK in "dsolow" "shock_SBT" {
	use datashocks_more
	local regressors = "`SHOCK'"
	forvalues lag = 1/`maxlag' {
		local regressors = "`regressors' L`lag'.`SHOCK'"
	}
	noi reg D.`VAR' `regressors' /* post97 */
	*pause 
	use irfs
	gen `VAR'_`SHOCK' = 0
	gen `VAR'_`SHOCK'_se = 0
	replace `VAR'_`SHOCK' = _b[`SHOCK'] if t==0
	replace `VAR'_`SHOCK'_se = _se[`SHOCK'] if t==0
	local sum = "0=`SHOCK'"
	forvalues lag = 1/`maxlag' {
		replace `VAR'_`SHOCK' = L.`VAR'_`SHOCK' + _b[L`lag'.`SHOCK'] if t==`lag'
		local sum = "`sum'+L`lag'.`SHOCK'"
		*disp "`sum'"
		save irfs, replace
		use datashocks_more
		test `sum'
		use irfs
		local tstat = sqrt(r(F))
		qui su `VAR'_`SHOCK' if t==`lag'
		replace `VAR'_`SHOCK'_se = abs(r(mean))/`tstat' if t==`lag'
	}
	replace `VAR'_`SHOCK' = L.`VAR'_`SHOCK' if `VAR'_`SHOCK'==0 & L.`VAR'_`SHOCK'~=0 & L.`VAR'_`SHOCK'~=.
	replace `VAR'_`SHOCK'_se = L.`VAR'_`SHOCK'_se if `VAR'_`SHOCK'_se==0 & L.`VAR'_`SHOCK'_se~=0 & L.`VAR'_`SHOCK'_se~=.
	gen `VAR'_`SHOCK'_lb = `VAR'_`SHOCK' - `VAR'_`SHOCK'_se
	gen `VAR'_`SHOCK'_ub = `VAR'_`SHOCK' + `VAR'_`SHOCK'_se
	save irfs, replace
}
}

foreach VAR of varlist *_dsolow {
	label var `VAR' "Response to Solow residual"
}
foreach VAR of varlist *_shock_SBT {
	label var `VAR' "Response to SBT shock"
}

pause on

foreach VAR in "sp" "prod" "hrs" {
	twoway (line `VAR'_shock_SBT t, lcolor(blue) lwidth(medthick) lpattern(solid)) /*
	*/ (line `VAR'_shock_SBT_lb t, lcolor(red) lwidth(medthick) lpattern(dash)) (line `VAR'_shock_SBT_ub t, lcolor(red) lwidth(medthick) lpattern(dash)) /*
	*/ if t>=0 & t<=10, yline(0, lwidth(thin) lcolor(black)) ytitle(" ", size(zero)) ylabel(-1(.2)1, angle(horizontal)) xtitle(" ") xtitle(, size(zero)) xlabel(0(2)10) /*
	*/ title("`VAR'", size(medium)) legend(off) name(`VAR'_shock_SBT, replace) graphregion(margin(zero) fcolor(none) ifcolor(none)) plotregion(margin(zero) fcolor(none) ifcolor(none))

	*pause
	twoway (line `VAR'_dsolow t, lcolor(blue) lwidth(medthick) lpattern(solid)) /*
	*/ (line `VAR'_dsolow_lb t, lcolor(red) lwidth(medthick) lpattern(dash)) (line `VAR'_dsolow_ub t, lcolor(red) lwidth(medthick) lpattern(dash)) /*
	*/ (line `VAR'_shock_SBT t, lcolor(green) lwidth(medthick) lpattern(dot)) /*
	*/ if t>=0 & t<=10, yline(0, lwidth(thin) lcolor(black)) ytitle(" ", size(zero)) ylabel(-1(.2)1, angle(horizontal)) xtitle(" ") xtitle(, size(zero)) xlabel(0(2)10) /*
	*/ title("`VAR'", size(medium)) legend(off) name(`VAR'_dsolow, replace) graphregion(margin(zero) fcolor(none) ifcolor(none)) plotregion(margin(zero) fcolor(none) ifcolor(none))
	*pause

	* Output data for IRFs for plotting in Matlab:
	outsheet `VAR'_shock_SBT `VAR'_shock_SBT_lb `VAR'_shock_SBT_ub using `VAR'_svar_`maxlag'lags.txt, nonames nolabel replace
	outsheet `VAR'_dsolow `VAR'_dsolow_lb `VAR'_dsolow_ub using `VAR'_solow_`maxlag'lags.txt, nonames nolabel replace

}
graph combine sp_shock_SBT sp_dsolow prod_shock_SBT prod_dsolow hrs_shock_SBT hrs_dsolow, colfirst rows(2) xsize(3) ysize(1.75) 
graph save solow_irfs_basicVAR_`maxlag'lags, asis replace
graph export solow_irfs_`maxlag'lags.eps, as(eps) replace logo(off) cmyk(off)
pause
*outsheet using irfs_solowres_basicVAR.csv, comma nolabel replace
*outsheet using irfs_solowres_basicVAR.txt, nonames nolabel replace

graph drop _all

foreach VAR in "sp" "wL" "wH" "rele" "hL" "hH" {
	twoway (line `VAR'_shock_SBT t, lcolor(blue) lwidth(medthick) lpattern(solid)) /*
	*/ (line `VAR'_shock_SBT_lb t, lcolor(red) lwidth(medthick) lpattern(dash)) (line `VAR'_shock_SBT_ub t, lcolor(red) lwidth(medthick) lpattern(dash)) /*
	*/ if t>=0 & t<=10, yline(0, lwidth(thin) lcolor(black)) ytitle(" ", size(zero)) ylabel(-1(.2)1, angle(horizontal)) xtitle(" ") xtitle(, size(zero)) xlabel(0(2)10) /*
	*/ title("`VAR'", size(medium)) legend(off) name(`VAR'_shock_SBT, replace) graphregion(margin(zero) fcolor(none) ifcolor(none)) plotregion(margin(zero) fcolor(none) ifcolor(none))

	*pause
	twoway (line `VAR'_dsolow t, lcolor(blue) lwidth(medthick) lpattern(solid)) /*
	*/ (line `VAR'_dsolow_lb t, lcolor(red) lwidth(medthick) lpattern(dash)) (line `VAR'_dsolow_ub t, lcolor(red) lwidth(medthick) lpattern(dash)) /*
	*/ (line `VAR'_shock_SBT t, lcolor(green) lwidth(medthick) lpattern(dot)) /*
	*/ if t>=0 & t<=10, yline(0, lwidth(thin) lcolor(black)) ytitle(" ", size(zero)) ylabel(-1(.2)1, angle(horizontal)) xtitle(" ") xtitle(, size(zero)) xlabel(0(2)10) /*
	*/ title("`VAR'", size(medium)) legend(off) name(`VAR'_dsolow, replace) graphregion(margin(zero) fcolor(none) ifcolor(none)) plotregion(margin(zero) fcolor(none) ifcolor(none))
	*pause

	* Output data for IRFs for plotting in Matlab:
	outsheet `VAR'_shock_SBT `VAR'_shock_SBT_lb `VAR'_shock_SBT_ub using `VAR'_svar_`maxlag'lags.txt, nonames nolabel replace
	outsheet `VAR'_dsolow `VAR'_dsolow_lb `VAR'_dsolow_ub using `VAR'_solow_`maxlag'lags.txt, nonames nolabel replace

}
graph combine sp_shock_SBT sp_dsolow wL_shock_SBT wL_dsolow wH_shock_SBT wH_dsolow rele_shock_SBT rele_dsolow hL_shock_SBT hL_dsolow hH_shock_SBT hH_dsolow, colfirst rows(2) xsize(3) ysize(1.75) 
graph save solow_irfs_morevars_`maxlag'lags, asis replace
graph export solow_irfs_`maxlag'lags_addvars.eps, as(eps) replace logo(off) cmyk(off)
pause
*outsheet using irfs_solowres_morevars.csv, comma nolabel replace
*outsheet using irfs_solowres_morevars.txt, nonames nolabel replace
save irfs_solowres, replace
outsheet using irfs_solowres.txt, nonames nolabel replace

graph drop _all

} /* end forvalues maxlag */



stop

* To get standard errors, do the Basu-Fernald-Kimball trick and estimate bivariate VARs
* Restrictions: shock is i.i.d., i.e. all lag coefficients zero (extended SR zero restriction)

use datashocks_more

set more off
pause on
irf set solow_irfs, replace

foreach SHOCK in "shock_SBT" "dsolow" {
foreach VAR in "sp" "prod" "hrs" "rels" {

	constraint define 1 [`SHOCK']LD.`VAR'=0
	constraint define 2 [`SHOCK']L2D.`VAR'=0
	constraint define 3 [`SHOCK']L3D.`VAR'=0
	constraint define 4 [`SHOCK']L4D.`VAR'=0

	constraint define 5 [`SHOCK']L.`SHOCK'=0
	constraint define 6 [`SHOCK']L2.`SHOCK'=0
	constraint define 7 [`SHOCK']L3.`SHOCK'=0
	constraint define 8 [`SHOCK']L4.`SHOCK'=0

	mat A = (.,0\.,.)
	mat B = (1,0\0,1)

	*noi var `SHOCK' D.`VAR', lags(1/2) const(1/2 5/6)
	*noi var D.`VAR', lags(1/2) exog(`SHOCK' L.`SHOCK' L2.`SHOCK')

	noi svar `SHOCK' D.`VAR', aeq(A) beq(B) lags(1/2) varconst(1/2 5/6)	

	irf create `SHOCK'_`VAR', step(20)

	*irf graph sirf, level(67) /* csirf ?? */
	*irf graph sirf, level(67) irf(`SHOCK'_`VAR') impulse(`SHOCK') response(D.`VAR')
	
	*pause

}
irf cgraph (`SHOCK'_sp `SHOCK' D.sp sirf) (`SHOCK'_prod `SHOCK' D.prod sirf) (`SHOCK'_hrs `SHOCK' D.hrs sirf) (`SHOCK'_rels `SHOCK' D.rels sirf), level(67)
pause
}













